100 REM  PHYSICAL EPHEMERIS OF MARS
105 REM 
106 DEFDBL A-Z
110 DT=57: REM  Value of ET-UT (sec)
111 REM    in 1990; change DT for
112 REM    other years as needed!!
115 P1=3.14159265#
120 P2=2*P1: RD=180/P1: DR=1/RD
125 GOSUB 330: REM  Get date and UT
130 D=(J-2451545!)+F
135 D=D+DT/86400!: REM  UT to ET
140 GOSUB 405: PRINT
145 READ A1,A2,D1,D2
150 T=T-1: REM  Cent from 2000.0
155 A1=A1+A2*T
160 D1=D1+D2*T
165 Q=A1: GOSUB 710: A1=Q*DR
170 Q=D1: GOSUB 710: D1=Q*DR
175 D9=COS(D5)*SIN(A1-A5)
180 N9=SIN(D1)*COS(D5)*COS(A1-A5)
185 N9=N9-COS(D1)*SIN(D5)
190 K=RD*ATN(N9/D9)
195 IF D9<0 THEN K=K+180
205 READ W0,W1,W2: GOSUB 304
215 REM
220 REM  Planetocentric decl. Earth
222 SD=-SIN(D1)*SIN(D5)
224 SD=SD-COS(D1)*COS(D5)*COS(A1-A5)
226 DE=RD*ATN(SD/SQR(1-SD*SD))
228 REM
230 REM  P.A. of Axis
232 SP=COS(D1)*SIN(A1-A5)
234 CP=SIN(D1)*COS(D5)
236 CP=CP-COS(D1)*SIN(D5)*COS(A1-A5)
238 PA=RD*ATN(SP/CP)
240 IF CP<0 THEN PA=PA+180
242 IF PA<0 THEN PA=PA+360
244 IF PA>=360 THEN PA=PA-360
248 REM
250 REM  Phase and magnitude
252 K4=((R4+RR)^2-R*R)/(4*R4*RR)
254 CI=(R4*R4+RR*RR-R*R)/(2*R4*RR)
256 I=RD*ATN(SQR(1-CI*CI)/CI)
258 M4=5*LOG(RR*R4)/LOG(10)
260 M4=-1.52+.016*I+M4
262 REM
270 A1$="Central meridian: ###.##"
272 A2$="P.A. of axis:     ###.##"
274 A3$="Decl. of Earth:   ###.##"
276 A4$="Angular diameter: ###.#"
278 A5$="Phase:            ###.###"
280 A6$="Magnitude:        ###.#"
282 PRINT USING A1$;C
284 PRINT USING A2$;PA
286 PRINT USING A3$;DE
288 PRINT USING A4$;9.359999/RR
290 PRINT USING A5$;K4
292 PRINT USING A6$;M4
300 END
304 REM   Rotations from epoch
305 W=W0+W1*(D-.0057755*RR)+W2*T
306 Q=W-K: GOSUB 710: C=Q
308 RETURN
309 REM      North Pole R.A.
310 DATA 317.681, +0.678
311 REM      North Pole Dec.
312 DATA 52.886,  +0.352
313 REM      Rotation
314 DATA 176.655,350.8919830,0.62
315 REM 
320 REM    Calendar --> JD
325 REM
330 INPUT "YR,MO,DAY  ";Y,M,D
335 INPUT "UT (H,M,S) ";H1,M1,S1
340 D=D+H1/24+M1/1440+S1/86400!
345 G=1: IF Y<1582 THEN G=0
350 D1=INT(D): F=D-D1-.5
355 J=-INT(7*(INT((M+9)/12)+Y)/4)
360 IF G=0 THEN 380
365 S=SGN(M-9): A=ABS(M-9)
370 J3=INT(Y+S*INT(A/7))
375 J3=-INT((INT(J3/100)+1)*3/4)
380 J=J+INT(275*M/9)+D1+G*J3
385 J=J+1721027!+2*G+367*Y
390 IF F>=0 THEN 400
395 F=F+1: J=J-1
400 RETURN
405 REM   Fundamental Arguments
410 REM
415 T=D/36525!+1
420 REM  T=Centuries from 1900.0
425 L0=.779072+.00273790931#*D
430 G0=.993126+.0027377785#*D
431 L4=.987353+.00145575328#*D
432 G4=.053856+.00145561327#*D
433 F4=.849694+.00145569465#*D
435 L5=.089608+.00023080893#*D
440 G5=.056531+.00023080893#*D
445 G6=.882987+9.294371E-05*D
450 L0=(L0-INT(L0))*P2
455 G0=(G0-INT(G0))*P2
456 L4=(L4-INT(L4))*P2
457 G4=(G4-INT(G4))*P2
458 F4=(F4-INT(F4))*P2
460 L5=(L5-INT(L5))*P2
465 G5=(G5-INT(G5))*P2
470 G6=(G6-INT(G6))*P2
475 REM      Sun
480 L=6910*SIN(G0)
485 L=L+72*SIN(2*G0)
490 L=L-17*T*SIN(G0)
495 L=L-7*COS(G0-G5)
500 R=1.00014-.01675*COS(G0)
505 R=R-.00014*COS(2*G0)
510 L=L0+L/206265!
515 REM      Mars
520 M=38451!*SIN(G4)
525 M=M+2238*SIN(2*G4)
530 M=M+181*SIN(3*G4)
535 M=M-52*SIN(2*F4)
540 M=M+37*T*SIN(G4)
545 M=M-22*COS(G4-2*G5)
550 M=M-19*SIN(G4-G5)
555 M=M+17*COS(G4-G5)
556 M=M+17*SIN(4*G4)
557 M=M-16*COS(2*G4-2*G5)
558 M=M+13*COS(G0-2*G4)
559 M=M-10*SIN(G4-2*F4)
560 M=M+7*COS(G0-G4)
561 M=M-7*COS(2*G0-3*G4)
564 B=6603*SIN(F4)
565 B=B+622*SIN(G4-F4)
570 B=B+615*SIN(G4+F4)
575 B=B+64*SIN(2*G4+F4)
580 R4=1.53031-.1417*COS(G4)
585 R4=R4-.0066*COS(2*G4)
590 R4=R4-.00047*COS(3*G4)
615 L4=L4+M/206265!: B4=B/206265!
620 N9=R4*COS(B4)*SIN(L4-L)
625 D9=R4*COS(B4)*COS(L4-L)+R
630 L1=ATN(N9/D9)
635 IF D9<0 THEN L1=L1+P1
640 LL=L1+L
645 V=N9*N9+D9*D9
650 RR=SQR(V+(R4*SIN(B4))^2)
655 S=R4*SIN(B4)/RR
660 BB=ATN(S/SQR(1-S*S))
665 E=(84428!-47*T)/206265!
670 N9=SIN(LL)*COS(E)-TAN(BB)*SIN(E)
675 D9=COS(LL): A5=ATN(N9/D9)
680 IF D9<0 THEN A5=A5+P1
685 S=SIN(BB)*COS(E)
690 S=S+COS(BB)*SIN(E)*SIN(LL)
695 D5=ATN(S/SQR(1-S*S))
700 RETURN
705 REM
710 REM    Normalize Degrees
715 Q=Q/360: Q=Q-INT(Q): Q=Q*360
720 RETURN
800 REM  This program computes 
810 REM  a number of interesting
811 REM  things about Mars' tele-
812 REM  scopic appearance at any
813 REM  date and time (UT).  Valid
814 REM  for several centuries, pro-
815 REM  vided that you reset
816 REM  the value of 'Delta-T' in
817 REM  line 110 to an appropriate
818 REM  value.  Written by Roger 
819 REM  Sinnott, this program is 
820 REM  explained more fully in 
821 REM  Sky & Telescope magazine
822 REM  for Sept. 1990, page 296.
